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ABSTRACT 

We characterize magnetically driven accretion at radii between 1 au and 100 au in 
protoplanetary discs, using a series of local non-ideal magnetohydrodynamic (MHD) 
simulations. The simulations assume a Minimum Mass Solar Nebula (MMSN) disc 
that is threaded by a net vertical magnetic field of specified strength. Confirming pre¬ 
vious results, we find that the Hall effect has only a modest impact on accretion at 
30 au, and essentially none at 100 au. At 1-10 au the Hall effect introduces a pro¬ 
nounced bi-modality in the accretion process, with vertical magnetic fields aligned 
to the disc rotation supporting a strong laminar Maxwell stress that is absent if the 
field is anti-aligned. In the anti-aligned case, we instead find evidence for bursts of 
turbulent stress at 5-10 au, which we tentatively identify with the non-axis 5 mimetric 
Hall-shear instability. The presence or absence of these bursts depends upon the de¬ 
tails of the adopted chemical model, which suggests that appreciable regions of actual 
protoplanetary discs might lie close to the borderline between laminar and turbulent 
behaviour. Given the number of important control parameters that have already been 
identified in MHD models, quantitative predictions for disc structure in terms of only 
radius and accretion rate appear to be difficult. Instead, we identify robust qualitative 
tests of magnetically driven accretion. These include the presence of turbulence in the 
outer disc, independent of the orientation of the vertical magnetic fields, and a Hall- 
mediated bi-modality in turbulent properties extending from the region of thermal 
ionization to 10 au. 

Key words: accretion, accretion discs - protoplanetary discs - MHD - instabilities - 
turbulence 


1 INTRODUCTION 


Gas in protoplanetary discs loses angular momentum and ac¬ 
cretes onto the central protostar at rates far greater than 
can be explained as a consequence of molecular viscosity 
( [Hartmann et al.|1998P . While myriad physical processes that 
might explain the observed rates continue to be debated (see 
the reviews by |Armitage|201l) [Turner et al.|20l4 l, the mag- 
netorotational instability (MRI; [Balbus & Hawleyj [19981 ) is 
widely considered to be the leading contender for providing 
the necessary enhanced angular-momentum transport. That 
this is true not just for protoplanetary discs, but also for a 
wide variety of accreting systems (e.g. compact objects), at¬ 
tests to the fact that the basic prerequisites for the instability 
are quite generic: a negative angular velocity gradient, a sub- 
thermal magnetic field, and enough free charges to provide 


sufficient coupling between the gas and the magnetic field. 
What makes protoplanetary discs particularly interesting from 
the perspective of the MRI is this final requirement; typical 
degrees of ionization in the deep interiors of such discs range 
from ~10“® to as low as and it is presently unclear 

how resilient the MRI is under these conditions. 

One of the first models of MRI-driven accretion in proto¬ 
planetary discs was constructed by |GamrnIe| l [ 1 996} , who pro¬ 
posed that Ohmic resistivity (due to weak coupling of charged 
species to the magnetic field) quenches the MRI in a ‘dead’ 
mid-plane region surrounded by ‘active layers’, which are ion¬ 
ized by cosmic rays. This layered-accretion scenario has been 
refined as additional non-thermal ionization sources (X-rays, 
UV radiation) and non-ideal magnetohydrodynamic (MHD) 
effects (ambipolar diffusion and the Hall effect) have been 
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mang, Terquem & Balbus[ 2002 

Salmeron & Wardle[[2003[ 

2005 2008 lllgner & Nelson[ 2006 Bai & Goodman[[2009| 

Bai & Stone||2011l [Wardle & Salmeron[12012[l. Several the- 

oretical studies, spanning both anal}rtical (e.g. [Blaes & Bal-[ 

bus|[ 19941 [Kunz & Balbus 20041 

Desch[ 20041 and numer- 

ical work (e.g. Hawley & Stone[ 

1998t [Bai & Stone[[2011[ 


|2013bt ISimon et al.j|2013b|al l, have shown that ambipolar 
diffusion plays a major role in determining the level of tur¬ 
bulence in protoplanetary discs. In the outer regions (radii 
R > 30 au), ambipolar diffusion acts to reduce the strength 
of MRI-driven turbulence close to the disc mid-plane; tur¬ 
bulence is stronger above this ‘ambipolar damping zone’ in 
a thin layer of gas strongly ionized by FUV photons ( [Perez-| 
[Becker & Chiang|2011|[slmon et al.|2013b|a| ). To produce ac¬ 
cretion rates consistent with observations, a large-scale mag¬ 
netic field threading the disc perpendicular to the disc plane is 
required ( [Simon et al.||2013b|a[ |. This net field enhances tur¬ 
bulence, despite the damping effect of ambipolar diffusion, 
and allows additional angular-momentum loss via a magnetic 
wind (akin to the [Blandford & Payne|198^ mechanism) . In the 
inner disc, both local (Bai & Stone 2013b I and global ( jGresselj 
jet aL|20f5] l simulations have shown that ambipolar diffusion 
can quench the MRI in the purported active layers of Gam- 
mie’s (1996) model. In those cases, a large-scale magnetic 
field is again necessary to drive accretion at the observation- 
ally inferred rates, this time solely via a magnetic wind in a 
manner akin to previous models of wind-driven accretion in 
protoplanetary discs (e.g. Wardle & Koenigl|1993t|Shu et al.j 


|1994t|Salmeron, Kdnigl & Wardle|2007l. 

The influence of the Hall effect has been considered in 
several studies (e.g. |Wardle||1999t [Balbus & Terquem||2001| 
jSano & Stone||2002a|b I, but only recently has it been pos¬ 
sible to numerically investigate the case in which the Hall 
term dominates over both the inductive term and the other 
non-ideal effects. The Hall effect operates at densities inter¬ 
mediate to those where Ohmic diffusion (high densities) and 
ambipolar diffusion (low densities) dominate (e.g. jWardlej 
|2007| |. Where dominant, the Hall effect leads to a qualita¬ 
tively different mode of angular-momentum transport ( jKunzj 
|& Lesur|2013t[Lesur, Kunz & Fromang||2014l|Bai||20T^ , and 

the nature of the associated stress depends on the orienta¬ 
tion of the large-scale magnetic field with respect to the an¬ 
gular momentum of the disc (as suggested by the linear anal¬ 
ysis of [Wardlejl 1 999| ). This is quantified by the sign of fl-B 
($2 is the angular velocity vector and B is the large scale 
magnetic field). For fl-B < 0, the Hall effect was found to 
reduce turbulent angular-momentum transport significantly; 
for fl-B > 0, the transport was found to be enhanced, occur¬ 
ring primarily through large-scale laminar stresses and mag¬ 
netic winds. 

In this paper, we further examine the role of magnetic 
fields and non-ideal MHD in driving accretion in protoplan¬ 
etary discs. While we include all three non-ideal effects - 
Ohmic dissipation, ambipolar diffusion, and the Hall effect - 
our main focus is on the influence of the latter. We first show 
that the Hall effect has minimal impact on the behaviour of 
the accretion flow in the outer disc. This confirms predictions 
based upon linear analyses as well as the recent numerical 
results of |Bai| ( [2015| ), who used a different numerical scheme 
that the one we employ herein. At smaller radii in discs with 
fl ■ B < 0, we identify a new regime of ‘bursty accretion’. We 


suggest that the physical mechanism for this bursty accretion 
is a Hall-mediated instability of whistler waves in a differen¬ 
tially rotating fluid. We also discuss the implications of the 
Hall effect for disc observables and the degree to which the 
Hall effect introduces a bi-modality in disc properties. 

The paper is organized as follows. In Section we de¬ 
scribe the equations of non-ideal MHD and the numerical al¬ 
gorithms we employ to solve them. We also detail our treat¬ 
ment of the ionization physics and chemistry and our choice 
of the simulation parameters used in our runs. Finally, we end 
this Section by defining the the relevant diagnostics computed 
from our simulations. In Sectionj^we present our results and 
analysis, focusing on the inner-disc and outer-disc regions 
separately. Section [^places our results within the context of 
current models for the structure and evolution of protoplane¬ 
tary discs and discusses potential observational implications. 
Our conclusions are summarized in Section!^ 


2 METHOD 

2.1 Shearing-box equations 


Our numerical simulations are carried out in the local 
shearing-box approximation ( [Hawley, Gammie & Balbusj 
[1995[ |, a useful framework for describing phenomena that 
vary on lengthscales much less than the large-scale properties 
of the disc. A small patch of the disc, co-orbiting with angular 
velocity f2 = Hodz at a fiducial radial location Ro, is repre¬ 
sented in Cartesian coordinates {x, y, z) with x = R — Ro and 
y = Ro<f> corresponding to the radial (R) and azimuthal {(f>) 
directions of a cylindrical coordinate system. Differential ro¬ 
tation is accounted for by including the Coriolis force and by 
imposing a background linear shear vo = —q^oxcy, where 


^ _ 1 dn{R) 

^ Ho d In i? 


R=Ro 


we use q = 3/2, appropriate for a Keplerian disc. 

In units such that the magnetic permeability is equal to 
unity, the equations of motions are the continuity equation, 

^ + V-{pv)^0, ( 1 ) 

the momentum equation, 

^ + V-{pvv -BB)+v(^P+ 

= 2pqD.oX — pVtoZ — 2Hodz X pv, (2) 


and the induction equation, 
dB 


dt 


VX (vxB) 


= - Vx 


VoJ + Vh 


JxB {JxB)xB 


B 


VA 


52 


(3) 


where p is the mass density, pv is the momentum density, B 
is the magnetic field, P is the gas pressure, and J — ^xB 
is the current density. We assume an isothermal equation of 
state, P = pCg, where Cs is the isothermal sound speed. From 
left to right, the source terms in Equation © correspond to 
radial tidal forces (gravity and centrifugal), vertical gravity, 
and the Coriolis force. The electromotive force (EMF) on the 
right-hand side of Equation ([^ includes contributions from 
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Ohmic diffusion (‘O’), the Hall effect (‘H’), and ambipolar dif¬ 
fusion (‘A’). The importance of these non-ideal terms is char¬ 
acterized by the corresponding diffusivity parameters ??o, ?)h, 
and T]A, respectively, whose magnitudes we specify in Section 

[Q] 


our model disc the column-density and temperature profiles 
of the minimum mass solar nebula (MMSN; |Hayashi||1981[ ), 
which are given respectively by 


E(i?) = 1700 


(— 

\ 1 au 


- 3/2 

g cm"^, 


(4) 


2.2 Numerical algorithm 

For our numerical calculations, we use a modified version of 
the PLUTO code l |Mignone et al. |2007l l. PLUTO integrates 
the shearing-box equations ([lJ-0 using a standard Go¬ 
dunov scheme with second-order-accurate spatial reconstruc¬ 
tion, a monotonised central flux limiter, and a second-order- 
accurate Runge-Kutta time integration. We use a modified 
HLL Riemann solver (see Appendix A of |Lesur, Kunz & Fro-| 
|mang|[2014l l to properly integrate the Hall effect in the in¬ 
duction equation (|^. The remaining non-ideal MHD terms. 
Ohmic diffusion and ambipolar diffusion, are computed using 
second-order-accurate finite differencing and added to the in¬ 
duction equation in a manner consistent with the method of 
constrained transport (CT; [Evans & Hawley|1988l l to preserve 
the solenoidality constraint V-B = 0 to machine precision. 
Ohmic and ambipolar diffusion are treated via a slightly modi¬ 
fied version of the resistivity algorithm of the public version of 
PLUTO to avoid repeated current calculations. The inclusion 
of the Hall effect requires slightly more care, as described in 
Appendix B of |Lesur, Kunz & Fromang| ( [2014p . For all of the 
non-ideal effects, we use the same implementation as those 
authors. 

The radial (x) boundary conditions are shearing periodic 
( [Hawley, Gammie & Balbus|1995P ; the azimuthal boundaries 
are periodic; and the vertical boundary conditions satisfy a 
modified-outflow condition, in which p is calculated based on 
hydrostatic equilibrium, Vx and Vy are extrapolated with zero 
vertical gradient, Vz is a standard outflow condition, and B is 
forced to be purely vertical. These boundary conditions have 
been shown to produce results similar to zero-current condi¬ 
tion on the boundary ( [Lesur, Kunz & Fromang|2014P . 

Shearing-box simulations of the MRI in vertically strati¬ 
fied discs generically exhibit significant mass outflow through 
the vertical boundaries of the computational domain (e.g. 
[Suzuki & Inutsul^|2009t [Simon et al.[||2013^ . Anticipating 
this, and with a desire to achieve a steady-state solution in 
our simulations, we continually replenish the disc mass by 
multiplying the density in each cell at each time step by a 
constant factor chosen such that the total mass within the 
computational domain is held constant. While this globally 
breaks conservation of momentum in our simulations - we 
do not decrease the velocity in each cell to compensate for 
the increased density - any local patch of a realistic accretion 
disc will, almost by definition, have mass inflow from larger 
radii. Our intention in replenishing the mass lost to vertical 
outflows is to numerically mimic this physical situation in a 
computational domain whose horizontal boundary conditions 
are shearing-periodic and thus preclude radial mass inflow. 


2.3 Simulation parameters 

Most of our ionization profiles and simulation parameters are 
taken from |Lesur, Kunz & Froman^ l |2014P , and we refer the 
reader to that work for more details. Briefly, we adopt for 


T{R) = 280 



(5) 


The density at the mid-plane, p{R), is determined by Equation 
l[^ and the relation E = s/^pH, where H{R) = Cs/U is the 
thermal pressure scale height of the disc at radius R. For a 
disc consisting of 80% H 2 and 20% He by number orbiting a 
solar-mass protostellar object, this gives 


p{R) = 1.4 X 10“® 



gem 


( 6 ) 


In this paper, we present simulations carried out at a series 
of radii, ranging from i?o = 1 au to Ro — 100 au and follow 
the dynamics up to a height 2 ; = ±6iTo, where Hq is the disc 
scale height at Ro. Previously published simulations by some 
of the authors ( [Lesur, Kunz & Fromang|2014P focused on the 
evolution of the MRI at Ro — 1, 5, and 10 au. We include 
some of these simulation data alongside our own. 

All of our model discs are initially in vertical hydrostatic 
equilibrium, 

P(»,y,2) = p(^?o)exp , (7) 


and are threaded by a uniform vertical magnetic field, B = 
Boez, whose strength and direction is determined by specify¬ 
ing the dimensionless free parameter 


_ Bo 2poCgo 

- jBoj^f^ 


( 8 ) 


where po is the initial density evaluated at the disc mid¬ 
plane. It is well established that the MRI, when subject to Hall 
EMFs, is sensitive to the relative orientation between the lo- 
cal magnetic field and the rotation axis (|Wardle[1999[[Balbus[ 


& Terquem 200l[ Sano & Stone[[2002a|b I ^We thus explore 

both magnetic polarities, with /3o = ±10'^, ±10^*, and ±10®. 
The linear growth of the MRI is seeded by perturbing the ra¬ 
dial velocity with a uniformly distributed random number be¬ 
tween ±0.1Ca. 

We quantify the diffusivities via several dimensionless pa¬ 
rameters that are independent of the magnetic-field strength 
(thus making them constant in time). For Ohmic resistivity, 
we define the magnetic Re}molds number 


Ru 


c,H 

VO 


(9) 


ambipolar diffusion is quantified by its Elsasser number, 

Am=^. (10) 

iloVA 


^ In fact, the Hall-MRI is also sensitive to the angle between the local 
magnetic field and the local vorticity l |Kunz|2008) ; for a Keplerian 
disc, the two of course coincide. 


© 2015 RAS, MNRAS 000,[lp6| 












































4 Simon et al. 



R (au) 


Figure 1. Representative ranges of Hm, Rhi 3tid Am at different radii 
in the mid-plane of our model disc, which is ionized by cosmic rays 
at a rate (^cr = 10“^^ exp(—S/96 g cm“^) jUmebayashi & Nakano] 
|1980| l and which contains dust grains (1% by mass) of varying radii 
(0.3 pm, 1 pm, 10 pm). Other chemistry parameters (e.g. grain ma¬ 
terial density, recombination rates, gas-grain collision rates) are taken 
from |Kunz & Mouschovias| l [2009) . 


Following |Lesur, Kunz & Fromang| ( [2014| ), we define the Hall 
Re 5 niolds number 

= ( 11 ) 

»7h «h 

where the Hall lengthscale £h is defined in situ. (In an ion- 
electron-neutral plasma, it is equivalent to the ion skin depth 
divided by the square root of the mass-weighted ionization 
fraction). 

In choosing the profiles of the magnetic diffusivities tjq, 
77h, and t]a, we are guided by the chemical model and gen¬ 
eralized Ohm’s law prescribed in §4 and appendix B of |Kunz| 
!& Mouschovias| ( [2009 1 and the setup used in the shearing- 
box simulations of Lesur, Kunz & Fromangj l |2014l l . In brief, 
for the simulations between i?o = 1 and 10 au, we adopt the 
same ionization structure as [Lesur, Kunz & Fromangj ( [2014P . 
At larger radii, we choose mid-plane values of i?m. Am, and 
Ru that are within the bounds calculated using the chemical 
model of jKunz & Mouschoviasj l |2009P with varying dust grain 
size; see Table [ij and Fig. [ij Note that in calculating our dif¬ 
fusivities, we have assumed that the mid-pla ne is ionized by 
cosmic rays ^cr = 10“^^ exp(—E/96 g cm“^) (Umebayashi & 


|Nakano|1980p . Recent results have brought into question the 
viability of cosmic rays as an ionization source, as there are 
concerns about geometrical effects significantly attenuating 
cosmic rays ( [Umebayashi & Nakano|2009P or the stellar w ind 
blocking them altogether ( Cleeves, Adams & Bergin[2013[ ?). 
Thus, we emphasize that our setup is just one model of many 
and that the cosmic ray flux used here should be taken as an 
upper limit. 

Away from the mid-plane, we adopt a simple approach 
to make contact with previous numerical setups (e.g. |Simon| 
[et al.jBOlBbT^ in which FUV photons are assumed to strongly 
ionize a thin layer above and below a diffusion-dominated 
mid-plane region. We set the diffusivities to be such that the 
corresponding Elsasser numbers (i.e. va/^oV> where va is the 
Alfven speed and rj is the corresponding diffusivity) are ini¬ 
tially constant (at their mid-plane values) for columns greater 


than 0.01 g cm“^ and assume a very high ionization fractio n 
(10“®) elsewhere (following Perez-Becker & Chiang 20111. 
This ionization depth corresponds to ~ 2Ho away from the 
disk mid-plane at these radii. In a Courant-Iimited integration 
in which diffusion is important, the time step is inversely pro¬ 
portional to the diffusivity and thus computations can become 
prohibitively expensive. In these regions of high ionization, 
the Hall and Ohmic Elsasser numbers are S> 1, but ambipolar 
diffusion, whose associated Elsasser number is proportional 
to gas density, can still be quite strong at large \z/H\. To min¬ 
imize the effect of strong ambipolar diffusion on reduci ng the 


time step, we cap rjA in these regions to values lOHo^o ■ Lesur, 


[Kunz & Fromangj ( [2014[ l examined the effect of changing this 
cap on the simulations at 1 au and found no significant differ¬ 
ences in their results; thus, we do not anticipate that the exact 
value of this cap will strongly influence the results here. In the 
simulations at 100 au, we do not include Ohmic resistivity at 
all, as it is expected to be negligible in these regions. 

Our simulations are conducted in units such that Cso = 
Tlo — Ho — po = For Ro = 1-10 au, the computational 


domain is chosen to have size Lj, x Ly x = 4 x 8 x 12, 
resolved by N^x NyX = 64 x 64 x 192 cells. At radii larger 
than 10 au, the domain size is L^x LyX = 8 x 16 x 12 and 
is resolved by x Ny x Nz = 128 x 256 x 192 cells. All of 
our simulations are listed in Table [l) following the convention 
of [Lesur, Kunz & Fromangj ( [2014| , the runs are labelled by 
their radial location in our model disc, the non-ideal MHD 
terms that are included, and by the strength and orientation 
of the magnetic field. For example, the shearing box in run 30- 
OHA-4n is placed at a disc radius Ro — 30 au, includes the 
Ohmic, ambipolar, and Hall terms, and has /3o = —10"^; the ‘n’ 
in 4n indicates that JT-B < 0 initially (as opposed to ‘p’ for 
Ct-B > 0). Two of our simulations used a mid-plane value of 
Am = 1 instead of Am = 10 as calculated using the methods 
described above. These runs are labelled with Ami’ appended 
to the simulation name. Finally, the simulations taken directly 
fromjLesur, Kunz & Froman^([2014[l are appended with ‘LKF’. 


2.4 Diagnostics 

We use several diagnostics to characterize the physics of ac¬ 
cretion in our shearing-box simulations. The first diagnostic is 
a volume average over the entire domain, 

(Q) = /// fi2) 

We also perform horizontal averages to plot quantities in 
(z, t) space-time diagrams; this average is denoted with an¬ 
gled brackets subscripted with xy. 

{Q).y = II dxdyQ. (13) 

Finally, we denote a time average with an over-bar, 

Q = ^ J dtQ, (14) 

where r is the time range over which we average. For each 
simulation, this averaging is chosen to begin after initial tran¬ 
sients have died out and to end when the simulation is termi¬ 
nated. 

Throughout much of this work, we examine the density- 
weighted i?</)-component of the stress tensor, pVxSvy — BxBy, 
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Table 1. Shearing-box simulations 


Label 

Radius 

(au) 

Rva 

Rh 

Am 

Po 


Q 

Qimid 

l-OHA-5p-LKF 

1 

2.9 

0.011 

0.24 

10 ® 

5.0 

X 

10 "® 

3.5 X 10 “® 

l-OHA-5n-LKF 

1 

2.9 

0.011 

0.24 

- 10 ® 

3.9 

X 

lO ""' 

- 9.0 X lO -®" 

l-OHA-3p-LKF 

1 

2.9 

0.011 

0.24 

10 ® 

3.1 

X 

10 "® 

2.2 

X 

10 “® 

5-OHA-3n 

5 

1400 

0.66 

2.5 

- 10 ® 

2.4 

X 

10 "® 

2.2 

X 

10 “® 

5-OHA-4n 

5 

1400 

0.66 

2.5 

- 10 '‘ 

5.0 

X 

10 "® 

8.4 

X 

10 “'® 

5-OHA-5p-LKF 

5 

1400 

0.66 

2.5 

10 ® 

2.0 

X 

10 "® 

2.1 

X 

10 “^ 

5-OHA-5n 

5 

1400 

0.66 

2.5 

- 10 ® 

3.3 

X 

10 "® 

2.3 

X 

10 “® 

5-OHA-5n-diff 

5 

1400 

0.66 

0.25 

- 10 ® 

6.3 

X 

10 "® 

3.1 

X 

10 “'® 

5-OHA-5n-axi 

5 

1400 

0.66 

2.5 

- 10 ® 

9.3 

X 

10 "® 

7.9 

X 

10 -" 

5-OA-5n 

5 

1400 

OO 

2.5 

- 10 ® 

9.2 

X 

10 "® 

1.0 

X 

10 “'® 

10-OHA-5p-LKF 

10 

10000 

1.8 

3.2 

10 ® 

1.0 

X 

10 "® 

9.3 

X 

10 “® 

10-OHA-5n 

10 

10000 

1.8 

3.2 

- 10 ® 

4.6 

X 

10 "® 

3.5 

X 

10 “® 

30-OHA-3p 

30 

43000 

6.5 

1 

10 ® 

4.0 

X 

10 "® 

9.8 

X 

10 “® 

30-OHA-3n 

30 

43000 

6.5 

1 

- 10 ® 

2.5 

X 

10 "® 

4.2 

X 

10 “'® 

30-OHA-4p 

30 

43000 

6.5 

1 

10 "^ 

9.3 

X 

10 "® 

1.2 

X 

10 “® 

30-OHA-4n 

30 

43000 

6.5 

1 

- 10 '‘ 

5.1 

X 

10 "® 

1.2 

X 

10 “® 

30-OHA-5p 

30 

43000 

6.5 

1 

10 ® 

2.5 

X 

10 "® 

7.4 

X 

10 “^ 

30-OHA-5n 

30 

43000 

6.5 

1 

- 10 ® 

1.3 

X 

10 "® 

8.0 

X 

10 “® 

100-HA-4p-Aml0 

100 

CX!) 

65 

10 

10 "^ 

5.2 

X 

10 "® 

2.3 

X 

10 “® 

100-HA-4n-Aml0 

100 

OO 

65 

10 

- 10 '‘ 

3.7 

X 

10 "® 

1.3 

X 

10 “® 

100-HA-4p-Aml 

100 

oo 

65 

1 

10 "^ 

1.3 

X 

10 "® 

9.0 

X 

10 “'® 

100-A-4p-Aml 

100 

OO 

OO 

1 

10 ‘‘ 

1.1 

X 

10 "® 

4.7 

X 

10 “^ 


which is responsible for the radial transport of angular mo¬ 
mentum. We cast this quantity in terms of the a parameter of 
IShakura & Syunyaev] ( [1973[ | by dividing by the gas pressure 
and performing the appropriate averages, 


(^PVx^'^y 


{pci) 


(15) 


The value of a for each simulation is given in Table 0. along 
with the same quantity integrated only over the disc mid¬ 
plane Qz\ < 0.5), Omid. 

while angular momentum can also be transported in 
the vertical direction by wind-launching mechanisms (e.g. 
[Blandford & Payne|1982[|Bai & Stone||2013a|bp , this angular- 
momentum transport is not well defined in a shearing box. As 
pointed out by several authors (e.g. Simon et al.|2013^|Bai &! 


|Stone||2013b[ |Lesur, Ferreira & Ogilvie||2013 1, the radial sym¬ 
metry of the shearing box can lead to magnetic-field struc¬ 
tures that launch outflows from the top and bottom of the 
box in opposite radial directions. Furthermore, the fast mag- 
netosonic point is not located within local domains ( |Lesur,| 
[Ferreira & Ogilvie[|2013l l. Therefore, we will not discuss this 
angular-momentum transport mechanism further. 



t (orbits) 


Figure 2. Volume-averaged total (Maxwell plus Reynolds) stress nor¬ 
malized by the volume-averaged gas pressure versus time. The black, 
solid line is run 5-OHA-5p (5 au, /3o = 10®) and the red, dashed line 
is 5-OHA-5n (5 au, /3o = —10®). The case with negative /3o shows 
significant temporal variations of up to an order of magnitude in am¬ 
plitude, whereas the positive Sq run shows very flat accretion stress 
after an initial transient phase. 


3 RESULTS 

3.1 Inner disc 

3.1.1 Simulations results 

In this Section, we consider radii 1-10 au, within which the 
Hall effect is the dominant non-ideal process (jWardle & Ngj 


1999t[Sano & Stone||2002a]|Wardle|2007l l. |Lesur, Kunz & Fro-j 


mangj (20141 showed that if fl-B > 0, the typical a. values 


are in the range 0 . 01 - 0 . 1 . Furthermore, the Maxwell stress 
is dominated by a laminar structure that exists at the largest 
horizontal scales of the box rather than by fluctuations that 
exist on smaller scales. [Lesur, Kunz & Fromangj l |2014l l also 
considered the case fl-S < 0 for i?o = 1 au (which we 
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Figure 3. Space-time diagram of the horizontally averaged Maxwell stress, normalized by the initial, mid-plane gas pressure for the two runs 
shown in Fig.[^ The top panel corresponds to 5-OHA-5p (5 au, /3o = 10®) and the bottom panel corresponds to 5-OHA-5n (5 au, Po = —10®). 
The Po < 0 case shows significant temporal variations in the Maxwell stress that occur near the disc mid-plane. The Po > 0 run shows a constant 
Maxwell stress in this same region. 


have included in our simulation suite) and found essentially 
no transport. 

We have expanded on these inner-disc simulations by 
carrying out simulations at 5 and 10 au with fl-B < 0. The 
evolution of the volume-averaged stress is shown in Fig. 
for two simulations at 5 au with Po = 10® (black, solid line) 
and Po = —10® (red, dashed line). In contrast to the stress 
measured in the il-i? < 0 simulations at 1 au, the stress mea¬ 
sured at 5 au undergoes significant temporal fluctuations for 
fl-B < 0, varying by up to an order of magnitude. Fig. 
shows the space-time evolution of the Maxwell stress for these 
same two runs. In both cases, the Maxwell stress is dominant 
near the disc mid-plane, but for fl-B > 0 the stress is al¬ 
most completely dominated by the largest horizontal scales 
in the box (i.e., — ky — 0), whereas the bursts of stress 

in the fl-B < 0 case have a non-negligible small-scale com¬ 
ponent. Given analytic expectations that an anti-aligned mag¬ 
netic field in an accretion disc subject to strong Hall diffusion 
should be MRI-stable (e.g. |Wardle|1999T l, such bursty turbu¬ 


lence with Omid ranging from up to is rather 

surprising. Evidently, even for fl-B < 0 in the Hall-dominated 
regime, it is possible to get significant amounts of turbulent 
activity and angular-momentum transport near the disc mid¬ 
plane. 

To better understand the conditions under which such 
variable accretion can occur, we performed several additional 
inner-disc simulations with fl-B < 0. We found that this vari¬ 
ability is also seen for /)o = —10® at 10 au, suggesting that the 
exact values of the diffusion coefficients do not matter so long 
as the Hall effect is the dominant term in the magnetic induc¬ 
tion equation and Am is larger than unity. We also checked 
the dependence of this bursty behaviour on the strength of the 
magnetic field. The equivalent 5 au run with Po = —10"^ (5- 
OHA-4n) shows similar variability, though the magnitude of 
fluctuations is less than an order of magnitude. Fig. shows 
the space-time evolution of the Maxwell stress for this run. 
The simulation spends a significant time in a state with very 
low stress near the mid-plane and highly variable stress in 
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Figure 4. Space-time diagram of the horizontally averaged Maxwell stress, normalized by the initial, mid-plane gas pressure for the run 5-OHA- 
4n. There appear to be two separate regimes of evolution. Until ~150 orbits, the Maxwell stress near the disc mid-plane is small, but in the upper 
atmosphere varies substantially. After 150 orbits, the variability in the active regions is reduced and there appear to be spasmodic increases in 
Maxwell stress within the mid-plane region. 


the upper atmosphere, followed by a more temporally fluctu¬ 
ating stress in the mid-plane starting after ~150 orbits. The 
simulation with /?o = —10^ at 5 au (run 5-OHA-3n) shows 
no Maxwell stress within the mid-plane region and significant 
large-scale stress for 1^1 >2. 

Finally, we recall that this large temporal variability was 
not seen in previous numerical studies of the inner-disc re¬ 
gion. For example, |Lesur, Kunz & Fromang| ( [2014P carried 
out fl-B < 0 simulations at 1 au and found the angular- 
momentum transport to be significantly reduced (a ~ 10“^ 
and no significant outflow). In their disc model at 1 au. Ohmic 
dissipation is as strong as the Hall effect, suggesting that the 
Hall effect must be dominant to trigger such bursty accre¬ 
tion. [B^ (j201^ explored Sl-B < 0 at 5 au, where the Hall 
effect does dominate over Ohmic dissipation, but used quasi- 
ID shearing boxes with negligible horizontal dimensions. This 
suggests that the turbulent activity we see in such a case re¬ 
quires either ^ 0, ky 0 or both. We test the conjecture 
that the burst activity requires non-axisymmetry (ky / 0) 
by carrying out a test simulation in which axisymmetry was 
enforced by setting Ny = 1 (run 5-OHA-5n-axi). This run 
showed no sign of bursty behaviour and the mid-plane region 
remained relatively quiescent. It is also worth noting that run 
5-OHA-3n, which does not exhibit such highly variable stress 
(for reasons that remain unclear), also does not exhibit signif¬ 
icant non-axisymmetric fluctuations. Follow-up work by [Bai| 
l [2015P did explore the case Cl-B < 0 at both 5 and 10 au in 
full non-axisymmetric 3D, but did not find the variability we 
see in our corresponding simulations. While [Bai| ( |2015P used 
a different chemistry calculation than the one employed here, 
resulting in the enhancement of all three non-ideal effects, we 
suspect that the absence of bursts in that work is specifically 
due to the enhancement of ambipolar diffusion, resulting in 
a relatively smaller value of Am compared to our work. To 
test this hypothesis, we carried out an identical simulation to 
5-OHA-5n but with the magnitude of tja amplified by a fac¬ 
tor of 10 everywhere; this run is labelled 5-OHA-5n-diff in 
Table We found no turbulent activity in the mid-plane for 


this test simulation. While it is unsurprising that a decrease in 
Am leads to less activity in the disc, it is important to point 
out that the value of Am that we have used in our fiducial 
simulations falls within the range of diffusion values based 
on our chemistry calculations (see Fig.[^. 

To summarize, we find that inner-disc regions of proto¬ 
planetary discs in which a net vertical magnetic field is anti¬ 
aligned with the rotation axis, the Hall effect dominates over 
magnetic induction and the other non-ideal processes, and 
neutrals collide at least a few times per orbit with free charges 
can support appreciable amounts of highly variable turbulent 
transport. What makes this conclusion all the more startling 
is that we find such disc regions to exhibit larger a than they 
would if the Hall effect were otherwise absent (see run 5-OA- 
5n in Table [^! To aid in our interpretation of this behaviour, 
we return to the linear theory. 

3.1.2 Origin of the turbulent bursts 

To provide some illumination on the observed bursty be¬ 
haviour seen in some of our simulations, we revisit the 
general linear analysis of the Hall effect in protoplanetary 
discs (cf. |Balbus & Terquem||20011 l. We begin by perturb¬ 
ing Equations liFW about the equilibrium state p = po, 
V — —gOoxtiy, and B = ByoCy -F B^oSz, the latter taken to 
be constant. Eulerian perturbations are denoted by a 5 and are 
taken to have space-time dependence cc5{t) exp(ik-x), where 
the time-dependent wavevector k — k{t) satisfies 

^ = qOokyez:. (16) 

To keep the analysis relatively simple, we shall ignore Ohmic 
dissipation and ambipolar diffusion and work in the incom¬ 
pressible limit at the disc mid-plane. The latter assumption 
precludes sound waves and (stable) buoyant oscillations. For 
notational convenience, we set Trrpo = 1, so that the Alfven 
velocity in the equilibrium state va,o = Bq. 

Under these assumptions. Equations Q-([^ written out 
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to linear order in 5 are, respectively, 

0 = k-Sv, 
dSv 


(17) 


dt 

dSB 

dt 


= —ikSn + 2QoSvyej: — (2 — q)Q.Q 5 vxky + iljaSB, 

= —qfloSBj;ey + iujaSv + inuJAkxSB, 


(18) 

(19) 


where we have introduced the compact notation loa = k-VAo 
for the Alfven frequency and 511 = 5P + Bo-5B for the 
perturbed total (gas + magnetic) pressure. Since dk^/dt = 
qQoky (Eq. |16| l. Equation ( |19| l together with Equation ( [17| l 
guarantee the divergence-free condition d{k-SB)/dt = 0. 
Note that the combination k-Bo is also time-independent. 

Motivated by the results described in Section Em we 
consider non-axis}mimetric fluctuations, ky / 0. To make 
analytical progress, we assume kx evolves slowly compared 
to the shear timescale which requires ky/kx <C 1. 

The coefficients in Equations l [18P and l [19P are then constant 
to leading order and we may consider solutions of the form 


S{t) = 5exp(7t). Using Equation (171 in the ^-component 
of Equation ( [ISP to determine the perturbed pressure, the x- 
and j/-components of the momentum equation ( [ISP can be 
then re-written in the following compact matrix form: 


lO^A 


7^ -I- k'^x “ 9)llo 


2flox'\ f^Bx 

7 ) \5By 


( 20 ) 


where = 2(2 — g) fig is the square of the epicyclic fre¬ 
quency and X = fcz /k^ is a geometrical factor accounting for 
the anisotropy of the inertial response. Using Equation ( [20P t o 
replace the perturbed velocity in the induction equation ( |19P , 
we obtain 


7 

gflo - fczWHX” 

, ^A 


kzVn 

7 

7 




2T!oX^ 


SB, 

SBz 


= 0, (21) 


= -K^X7^ + ( Ho - 


where 

7 ^ = 7 ^ -I- -1- a;i - kzVnqO.0. 

The corresponding eigenvectors are 

e 


where «h = is a scale-dependent Hall velocity. 

Setting the determinant of the 2x2 matrix in equation 
to zero results in the dispersion relation, which may be 
conveniently written as 

7 ^ -I- i7fcDH -I- ^ 7 ^ - i7fcDH -I- - kzVuqQo (^7^ + 

( 22 ) 

(23) 

(24a) 
(24b) 
(24c) 


5Vy = 


T 

— I 

LOA 7 ^ 


(2 - g)no7^ - Loi ^ 


SBx 


^By 

loa 


i^x 

Ly2 


7 + kzVH - 

kzVH 


2Ho - 

kzVH 

X 


kzVu 

X 


7 2 flo - 


where ^x = Svxlx is the radial displacement of a fluid ele¬ 
ment. We examine different limiting cases of the dispersion 
relation \22\ and the eigenvectors (eq.|24P, not only to make 


contact with the existing literature but also to elucidate the 
various physical processes at play. 

Let us first take the limit fig —>■ 0 (i.e. no rotation, no 
shear), in which case Equation \22\ becomes 


7 ^ =p i-ykvn -I- tji = 0 . 


(25) 


The positive-frequency solutions (w = —i 7 > 0) are given by 
Mh 


LO 

— = T 
loa 


+ 1 /1 + 


kin 

~Y~ 


(26) 


Eor small wavenumbers (klu <C 1), both of the above solu¬ 
tions reduce to Alfven waves; at large wavenumbers (fc^H S> 
1), right-handed waves (plus sign) go over to the high- 
frequency whistler-wave branch (a; « feun), whereas large-fc 
left-handed ion-cyclotron waves (minus sign) are cut off at a 
frequency 


OJh = 


^A 


VA,Q 

& 


k-Br 


kBn 


(27) 


With solid-body rotation included, the dispersion relation 
reduces neatly to 


(28) 


7'= =p i7A:t;H ( 1 - ) +oj\ + 2 Q.okzVn = 0, 

kzVn ) 


and the solutions 


become 


lua 2 y \ 2 J ujYi k 


(29) 


where £h = fe(l — 2Qox/kzVn) is the Hall lengthscale mod¬ 
ified by rotation - its presence signals the interplay between 
the dynamical epicycles in the bulk neutral fluid and the ‘mag¬ 
netic epicycles’ (i.e. circular polarization) induced by the Hall 
effect. Depending upon the orientation of the mean magnetic 
field, this circular polarization can increase or decrease the ef¬ 
fective Coriolis force; e.g. for fl-B > 0, the dynamical epicy¬ 
cle is slowed by Hall currents. In this case, the return magnetic 
tension force is effectively increased, a consequence of the last 
term in the square root in Equation \29) . (See §3 of [Balbus~&| 
|Terquem|2001| for an alternate discussion of this case.) 

Next, we restore the differential rotation and take the 
ideal-MHD limit (wh —>■ 0). Equation (221 then reverts to the 
standard MRI dispersion relation of|Balbus & Hawl^([1991|l, 




-«^x I 


( 7 ^ +c^i) +4HgXt^i, 


(30) 


for which unstable solutions exist if lo\ ~ 2gHgX < 0. Non- 
axisymmetric shearing waves are amplified strongest by the 
MRI when x ~ 1. As these waves sweep from leading to 
strongly trailing, x 0 and the final (destabilizing) term be¬ 
comes smaller and smaller. Eventually, the non-axisymmetric 
MRI is stabilized when x < Lo\/2qQ^. Hence, the ideal non- 
axisymmetric MRI is a transient instability (iBalbus & Hawley] 

Finally, including both the Hall effect and differential ro¬ 
tation leads to two additional instabilities, which are phys¬ 
ically distinct from the MRI. The first is the Hall-shear in¬ 
stability (HSI), which may be readily obta ined by letting 
(loa/XIo)^ ~ Rii/y/P <C 1 in Equation (211, thereby sever¬ 


ing the dynamical link between the induction and momen- 
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turn equations ( [Rudiger & Kitchatinov||2005j jPandey & War-| 
|dle|2012] l. The dispersion relation ( |22| l then reads 


(y'^ + tPxj ( 7 ^ + - k^vnqt^o'j 


= 0 . 


(31) 


The first factor in parentheses represents epicyclic motion 
in the bulk neutral fluid, an inherently stabilizing influence 
in a Keplerian disc. The second factor in parentheses repre¬ 
sents whistler waves, which are physically decoupled from 
the epicycles and may be driven unstable if |gf2o/<^H| > 1, 
i.e. if the time for an Alfven wave to travel a distance £h is 
longer than the time for a magnetic perturbation to grow by 
the Keplerian shear. In this case. Equation ( |24| l gives 
- 1/2 

~ 0(1). (32) 


SBtj 


qQo 


kzVn 


This instability, and its non-axisymmetric guise, was also 
found by jKunzj l |2008P when specializing the linear analysis to 
planar (i.e. non-rotating) shear flows. It is important to note 
that the stabilizing term, k'^v'^, for a shearing wave eventually 
dominates over the destabilizing term, kzVnqtlo- Hence, the 
non-axisymmetric HSI is also a transient instability. 

The second instability may be extracted by setting 
kz/kx = £ <€. 1, so that x ~ This asymptotic ordering pre¬ 
cludes both the MRI and the HSI from entering the analysis, 
as it puts the focus on strongly leading/trailing disturbances. 
Employing the Ansatz y ~ 0(e), we find that the dispersion 


relation 1221 at 0(1) may be factored as 

-I- \y\ + 2flofczWH] \^‘a + (2 - q)Q.okzVii^ = 0. (33) 
For q > 0, the perturbations are unstable when 


ujI 


< kzVu < — 


uj\ 


(34) 


(2 - q)Q.o " ' 2Ho 

which, if satisfied for one kx, is always satisfied. For this rea¬ 
son, non-axisymmetric shearing waves are linearly unstable 
without boundj^In this case. Equation 1241 gives 

1/2 

kzVn + 


SB: 


yjx 

' 2tlo 


CJa 


(2 - g)no 


kzVn -h 


2flo 


0(e). 


(35) 


In the limit Ho —> 0, this dispersion relati on r ecovers the ion- 
cyclotron wave, 7 ^ ~ —uiXlk^v^ (see Eq. 271; thus, the insta¬ 


bility arises from the interaction between ion-cyclotron waves 
and epicyclic motion in a differentially rotating disc. The ax- 
isymmetric version of this particular instability is referred to 
by [Pandey & Wardlej ( |2012| l as the ‘diffusive MRI’, because 
the growing azimuthal field is generated not by the differen¬ 
tial shear of the radial field (as for the standard MRI and the 
HSI) but rather by Hall currents. 


3.1.3 Which instability might be responsible for the observed 
turbulent bursts? 

The linear theory elucidated, we now return to those simula¬ 
tions exhibiting bursts of turbulent activity and analyze them 

^ Recall that we have ignored the effects of Ohmic dissipation and 
ambipolar diffusion in our linear analysis, which would attenuate the 
growth of strongly leading/trailing shearing waves. 


with the aim of determining whether any of the instabilities 
described in the preceding section might be associated with 
the bursts. Our starting point is Fig.|^ which provides pseudo¬ 
color plots of the X-, y-, and ^-components of the magnetic- 
field fluctuations SB — B — {B)xy (where the horizontal 
average is performed at the disc mid-plane, 2 = 0) during 
a strong burst (orbit number 157) in run 5-OHA-5n; at this 
time and vertical location, (Bx)xy = — 0 . 021 , (By)xy = 0 . 26 , 
and (Bz)xy — — 0 . 0045 . The top row shows these fluctua¬ 
tions in the x-y plane at 2 = 0; non-axisymmetry is readily 
apparent, with ky/kx > 0 . The bottom row shows the same 
fluctuations in the y-z plane at 2 : = —1; the dominant modes 
appear to be tilted in this plane, with ky/kz > 0 . This mode 
geometry implies kzVnqtlo — kzkyiuByqtlo > 0, which ex¬ 
onerates the non-axisymmetric diffusive MRI (see Eq. |34| ) . 
To test whether the non-axisymmetric HSI, which requires 
kzVnqtlo > 0, is responsible for the bursty behaviour, we 
ask whether (wa/Hq)^ ^ 1 for the observed fluctuations. For 
that, we perform two-dimensional Fourier transforms of the 
panels In Fig.[^ find the dominant (kx, ky,kz) mode, and use 
the aforementioned mean-field values to compute k-{B)xy. 
This procedure gives (kx,ky,kz) — tv x ( 1 , 1 / 4 , 5 / 6 ) and so 
(cja/Hq)^ = 0 . 016 , which is indeed much smaller than unity. 
We also check whether these fluctuations satisfy the HSI po¬ 
larization (Eq. |32p . Using kzVn = 0.50 and x = 0 . 40 , as 
measured for the largest-amplitude fluctuations, in Equation 
l |32p , we can predict SBx/SBy = —0.70. A comparison of 
this value with the average value measured in Fig. namely 
SBxjSBy = —0.70, strongly supports our conjecture. 

While we have not definitively proven that the non- 
axisymmetric HSI is responsible for driving the large- 
amplitude variability seen in some of our simulations, we 
have shown that the fluctuations observed in these simula¬ 
tions that are associated with this variability have attributes 
consistent with those expected from this instability. In any 
case, the important point remains. Under certain conditions 
pertinent to the inner regions of protoplanetary discs, the 
Hall effect, acting on a differentially rotating disc that is 
threaded by a weak magnetic field antl-aligned with the ro¬ 
tation axis, can lead to bursts of turbulence and enhanced 
angular-momentum transport. These are the first disc simula¬ 
tions to demonstrate this. 

Finally, as has been pointed out by several authors (e.g. 
Simon et al. |2013a[ [Bai & Stone||2013b[ [Lesur, Ferreira &| 
Ogilvie|2013 1, outflows can be launched from the disk at the 
top and bottom of the shearing box, but in opposite radial 
directions. In this (unphysical) geometry, the toroidal field at 
any given time does not reverse sign along the 2 dimension, 
as it would If the outflows had a physically realistic geome¬ 
try. When present, the current sheet associated with this sign 
flip could lead to additional dissipation that would quench 
the non-axisymmetric HSI modes and prevent the bursts from 
occurring in the first place. We have examined the toroidal 
field geometry for the runs that demonstrate bursts and find 
that while two of the simulations have no toroidal field re¬ 
versal (5-OHA-4n and 5-OHA-5n), the remaining simulation 
(10-OHA-5n) does have such a reversal and still undergoes 
this variable accretion. Thus, we do not expect that the large 
scale magnetic field geometry will affect the existence of these 
bursts. 
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Figure 5. Magnetic-field fluctuations in (top row) the x-p plane at 2 = 0 and (bottom row) the y-z plane at a; = — 1 at orbit number 157 in 
run 5-OHA-5n, during the largest-amplitude burst of Maxwell stress. From left to right, the panels show the x-, y-, and 2 -components of the 
fluctuating magnetic field SB = B — {B{x,y,z = 0))a:y', {Bx)xy = —0.021, (By)xy = 0.26, and {Bz)xy = —0.0045. The dashed line in the 
top (bottom) row indicates the planar cut shown in the bottom (top) row. Th e fluctu ations exhibit non-axisymmetry with > 0 and 

(oj/O. q)'^ <g 1, implicating the non-axisymmetric HSI as the source of stress ( |3.1.3| l. These characteristics are generic to all of the bursts that 
we observed in this run. 


3.2 Outer Disc 


In this Section, we focus on radii larger than 10 au. Fig. 
shows the space-time evolution of the horizontally averaged 
xy Maxwell stress for runs 30-OHA-5p and 30-OHA-5n (/3o = 
± 10®). There is very little to no Maxwell stress within the mid¬ 
plane region of 30-OHA-5n (i.e. below the FUV-ionized sur¬ 
face layers, which corresponds to 1^1 < 2Ho), whereas there 
is significant stress in this region for 30-OHA-5p. This result 
is in agreement with recent numerical work by ( [Bai 20151 
and with the expectation, borne out of linear theory ( Wardle| 
|1999t [Balbus & Terquem||20011 l, for the Hall effect to make 
the disc more (less) unstable when fl-B > 0 (f2-B < 0). 
Indeed, there is nearly two orders of magnitude difference 
in Omid between the /3o = 10® and /3o = —10® cases. Other 
initial magnetic-field strengths (see Table show similar re¬ 
sults, though the differences between Omid for the two field 


polarities are less dramatic. The total volume-integrated a, on 
the other hand, decreases by just a factor of order unity when 
the polarity of the net vertical field is reversed. In summary, 
the Hall effect affects the mid-plane djmamics significantly at 
30 au, but does not have a large impact on the total integrated 
stress suggesting that the total volume-integrated stress is 
dominated by the surface contribution, consistent with the 
strong stress observed at high \z\ in Fig.|^ 

Fig.[^shows the space-time evolution of the horizontally 
averaged Maxwell stress from the simulations at 100 au with 
Am = 1 at the mid-plane with (100-HA-4p-Aml) and with¬ 
out (100-A-4p-Aml) the Hall effect. Comparing Omid for these 
two runs, we see that the mid-plane stress differs by less than 
a factor of two. We also examined the effect of increasing the 
mid-plane value of Am from 1 to 10 (runs 100-HA-4p-Aml0 
and 100-HA-4n-Aml0). Fig. gives a space-time compari- 
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Figure 6. Space-time diagram of the horizontally averaged Maxwell stress, normalized by the initial, mid-plane gas pressure for runs 30-OHA-5p 
(30 au, /So = 10®; top) and 30-OHA-5n (30 an, /So = —10®; bottom). While both runs show significant Maxwell stress for \z\ > 3, the /So < 0 
simulation has very little stress in the mid-plane region while the /So > 0 case exhibits a significant Maxwell stress in this region. 


son of the Maxwell stress in these runs, which have different 
magnetic-field polarities. From this plot, and by comparing 
the values of ctmid between 100-HA-4p-Aml0 and 100-HA- 
4n-AmlO, we find that the stresses in the mid-plane region are 
mostly independent of the field orientation. Ambipolar diffu¬ 
sion, not the Hall effect, is largely responsible for the stress 
at these large radii (consistent with the recent work of |Bai| 

|20T^ . 

3.3 Summary of results 

To summarize Section we have plotted a versus R for 
(nearly) all of our simulations in Fig. The strength of an¬ 
gular momentum transport strongly depends on the magni¬ 
tude and orientation of the magnetic field with respect to the 
disc angular momentum. The dependence on the orientation 
weakens significantly towards large radii and in simulations 
that undergo ‘bursts’ of turbulence (at R — 5-10 au ). 


4 IMPLICATIONS 

Our results have implications both for theoretical modeling 
and for future observations. Here, we discuss these impli¬ 
cations and construct an overall scenario for mass accretion 
in magnetized protoplanetary discs based on our results and 
those of other recent numerical efforts. 

The most robust inference from current simulations is 
that, in magnetized protoplanetary discs, a vertical magnetic 


field threading the disc is necessary to generate rates of 
angular-momentum transport consistent with observational 
constraints. This result applies to both the inner and outer 
discs ( |Bai & Stone||2013bi |Bai||2013] [Simon et al.||2013b|a| ), 
and immediately implies that the dependence of the Hall ef¬ 
fect on the magnetic-field direction (e.g. |Wardle|1999] ) might 
have an observable effect. 

The nature of the magnetic stresses responsible for 
angular-momentum transport varies with radius via changes 
in the net magnetic-field strength, the chemical properties of 
the disc, and the degree and nature of the disc ionization. 
All three of the mechanisms that are possible in principle - 
turbulence, laminar stresses internal to the disc, and large- 
scale wind torques - are found to be important in simulations. 
In the inner disc (generally R < 10-30 au), [Bai| ( |2015| l and 
[Lesur, Kunz & Fromangj l |2014l l have emphasized the domi¬ 
nant role of laminar stresses. The MRI can be suppressed at 
these radii, with the internal laminar MHD stresses arising 
from the interaction of Keplerian shear and the Hall effect 
(HSI). Such laminar stresses should be interpreted with cau¬ 
tion, as they have so far been seen only in local shearing- 
box simulations, which are not designed to capture features 
with horizontal correlation lengths larger than the disc verti¬ 
cal scale height. Those generated by the HSI in our simula¬ 
tions have effectively infinite horizontal correlation lengths, 
and so the simulations only indicate that the disc is suscep¬ 
tible to generating large-scale global structures. The nature 
of these structures, and their modeling within the confines of 
a-disc theory, remain largely speculative. 
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Figure 7. Space-time diagram of the horizontally averaged Maxwell stress (normalized by the initial, mid-plane gas pressure) for runs lOO-OHA- 
4p-Aml (top) and 100-OA-4p-Aml (bottom). The neglect of the Hall effect in run 100-OA-4p-Aml does not affect the evolution of the Maxwell 
stress. 
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Figure 8. Space-time diagram of the horizontally averaged Maxwell stress (normalized by the initial, mid-plane gas pressure) for runs lOO-OHA- 
4p (100 an, jSo = 10"*; top) and 100-OHA-4n (100 an, /3o = —10^; bottom). Both runs show significant stress near the disc mid-plane, which is 
independent of the orientation of the initial vertical magnetic field. 
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In the outer disc (R > 30 au), the stresses responsible 
for accretion can either be turbulent or laminar, depending 
upon the strength of the magnetic field and the column depth 
to which FUV photons ionize the gas ( [Simon et al.[|2013b|aP . 
In the case of turbulent accretion, the turbulence is ‘layered’, 
similar to the classic dead-zone model of |GammIe| l |1996| l, but 
with strong ambipolar diffusion near the mid-plane signifi¬ 
cantly reducing the efficacy of magnetorotational turbulence. 
Because this turbulence is not fully quenched in this region, it 
has been referred to as the ‘ambipolar damping zone’. 

The results presented here have added to and expanded 
upon this emerging paradigm for disc accretion in magne¬ 
tized protolanetary discs. One general consideration, which 
we have emphasized here, is that our ability to accurately 
model MHD stresses depends upon a knowledge of the rela¬ 
tive importance of different non-ideal effects as a function of 
radius and height within the disc. The general level of agree¬ 
ment between our results and those of [Baij ( |2015| ) suggests 
that the overall scenario outlined above is robust to algorith¬ 
mic differences in the treatment of the Hall effect. However, 
within the range of diffusion coefficients informed by chem¬ 
ical models, there are regions of parameter space in which 
(i) turbulence can be induced in the inner-disc region if the 
large-scale magnetic field is anti-aligned with the angular mo¬ 
mentum vector of the disc and (ii) turbulence can be en¬ 
hanced near the mid-plane in the outer disc due to lower 
ambipolar diffusion. Our results thus suggest that discs are 
not far from the boundary between laminar and turbulent be¬ 
haviour, and that generally similar discs could fall either way 
depending on the exact details of the ionization and recom¬ 
bination processes at work. While improving our modeling 
of ionization/recombination physics through, e.g., improved 
chemistry calculations and dust modeling will help us con¬ 
verge on a more robust set of predictions for the behaviour 
of accreting gas, perhaps the most promising avenue forward 
is to develop stringent observational tests of our theoretical 
models and use those tests to constrain the behaviour of ob¬ 
served disc systems. 

Recently, a promising avenue for testing accretion the¬ 
ories has emerged, namely to predict and then observation- 
ally constrain signatures of turbulence in protoplanetary discs 
IjCarr, Tokunaga & Najita||2004| [Hughes et aL||201l[ [Simon, 
Armitage & Beckwith [2011 Guilloteau et al. [2012 Forgan, 


Armitage & Simon| 2012 de Gregorio-Monsalvo et ai.[[2013| 

Shi & Chiang[[2014[ [Simon et al.[[2015[ |. As illustrated in 
Fig. 10 our results suggest qualitatively that the magnitude 
and nature of angular momentum transport are bi-modal, de¬ 
pending upon the orientation of the net vertical magnetic field 
with respect to the rotation axis. If fl-B > 0, the stress be¬ 
tween i? ~ 1 au and i? ~ 10 atj^is large and mostly laminar; 
if fi-B < 0, the stress is on average smaller and may ex¬ 
hibit bursty behaviour. The Hall effect, which is responsible 
for this behaviour, becomes less important farther out in the 
disc, as we have shown here, consistent with both anal}Aic 
( [Salmeron & Wardle|2005| ||2008[ l and numerical ( [Bai|2015] l 
expectations. 

[Simon et al.[ ( [2015[ | carried out numerical simulations 
focused on the outer regions of protoplanetary discs, which 


^ The exact radial locations that encompass these different regimes 
will depend on the surface density model adopted. 



Figure 9. a versus radius as measured from our shearing-box simu¬ 
lations. The colors denote the magnetic-field strength while the sym¬ 
bols denote polarity, with squares (asterisks) corresponding to /3o > 0 
(< 0). Multiple values of a at a given radius result from different 
simulations performed at the same radius but with different diffusion 
values (see Table[^. All runs except for l-OHA-5n-axi, 5-OA-5n, and 
100-A-4p-Aml are included in this plot. At large radii, the magnetic 
polarity has minimal effect on the value of a, unlike at small radii. 
The presence of turbulence ‘bursts’ (circled points) reduces the effect 
of the magnetic polarity on the amplitude of the stress. 


included Ohmic and ambipolar diffusion, in order to pre¬ 
dict signatures of turbulence that could be constrained via 
molecular-line observations. Our results suggest that the pre¬ 
dictions made by [Simon et al.[ ( [2015| are not substantially af¬ 
fected by their neglect of the Hall effect at radii R ~ 100 au. 
Sub-mm observations are therefore predicted to show that 
the large-scale properties of discs are consistent with a single 
population. However, observations sensitive to scales smaller 
than 30 au could start to reveal the bi-modality produced in 
Hall-dominated accretion. Observations in the infrared ( [Carr,[ 
[Tokunaga & Najita[2004[ l are a potentially better probe of tur¬ 
bulence in the inner disc regions and could thus be a very im¬ 
portant tool for testing the predicted bi-modality of turbulent 
behaviour. 


5 CONCLUSIONS 

We have presented results from a series of local stratified 
shearing-box simulations of magnetorotationally driven pro¬ 
toplanetary disc accretion. The simulations assume a disc 
model based on the minimum mass solar nebula l [Hayashi| 
[1981[ l, in most cases threaded by a weak flux of vertical mag¬ 
netic field. Our main conclusions are: 

(i) The mid-plane of protoplanetary discs is not necessar¬ 
ily ‘dead’ and may have either turbulent or laminar Maxwell 
stresses that drive accretion at an enhanced rate (a ~ 0.001- 
0.1). 

(ii) In the inner disc (7?o < 30 au), the nature of the mid¬ 
plane angular-momentum transport depends upon the orien¬ 
tation of the large-scale vertical magnetic field with respect 
to the disc rotation axis, Q,-B. If the field and rotation vec¬ 
tors are aligned (i.e. f2-B > 0), angular momentum is trans¬ 
ported through a laminar Maxwell stress, in agreement with 
previous simulations by [Lesur, Kunz & Fromang[ l|2014[l and 
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Figure 10. Illustration of the predicted structure of poorly ionized protoplanetary discs threaded by weak (/3o ~ ±10'*“®) vertical magnetic fields. 
The structure of the inner, thermally ionized disc, and the outer disc where ambipolar diffusion is dominant, is dependent upon the strength of 
the net field but is independent of its sign. Between ~1 an and 10 an, however, the Hall effect is dominant and the disc structure depends on the 
sign of fl-B. 


|Bai| ( [2015| ). If the field is anti-aligned with the rotation axis 
(i.e. fl-B < 0), then the mid-plane either may remain qui¬ 
escent or may undergo ‘bursts’ of turbulent activity, depend¬ 
ing upon the precise values of the diffusivities and magnetic- 
field strength in that region. We presented evidence that these 
bursts of turbulence are driven by the non-axisymmetric ver¬ 
sion of the HSI, which results from the interaction of Kep- 
lerian shear and non-axisymmetric whistler waves on a pre¬ 
dominantly toroidal magnetic field. 

(hi) At intermediate radii, Rq ~ 30 au, the sign of 
Ct-B determines whether or not there is angular-momentum 
transport within the mid-plane region. However, the height- 
integrated stress is not strongly dependent upon the 
magnetic-field orientation - a result of the FUV ionization 
layer penetrating closer to the mid-plane at large radii. 

(iv) At large radii, Ro > 100 au, the orientation of the mag¬ 
netic field has little influence on the turbulent activity in the 
disc mid-plane and on the overall stress. This corroborates 
the conclusions of |Bai| ( [2015P , who used a different numeri¬ 
cal scheme for solving the Hall-MHD induction equation. 

A long standing goal of accretion-disc theory is to predict 
the structure and evolution of the disc from a first-principles 
understanding of the stresses that drive accretion. Our work 
does not accomplish that goal, first because we assume a par¬ 
ticular radial disc structure and chemistry model in calculat¬ 
ing the non-ideal terms and, secondly, because we do not ad¬ 
dress the transport due to winds, for which global simula¬ 
tions are needed. However, local simulations are now rela¬ 


tively computationally inexpensive, and it is possible to cal¬ 
culate the internal stresses as a function of the known control 
parameters (e.g. radius, surface density, external magnetic- 
field strength, disc chemistry) by carrying out a large num¬ 
ber of simulations that sweep through parameter space. The 
simulations carried out in this work represent a small, yet in¬ 
formative, region of this parameter space. While much of this 
parameter space remains to be probed via numerical simula¬ 
tions, our results suggest that any universal prediction to be 
made for quantities such as E(M, R) will be quite difficult to 
formulate, at least if MHD stresses dominate the evolution. 
Instead, quite minor changes to the adopted chemistry model 
can make a qualitative difference to the angular-momentum 
transport in some regions of the disc. Moreover, the net mag¬ 
netic field is a difficult-to-measure quantity, which will likely 
differ from disc to disc and which will evolve over time in a 
complex way dependent upon the radial structure of the disc 
turbulence. An MHD analog of an a-model for protoplane¬ 
tary discs, which predicts the structure given a small number 
of parameters that could be observationally constrained, does 
not seem realistically attainable at this time. 

Adopting a more optimistic tone, there are excellent 
prospects for testing whether the MHD processes that we have 
discussed are actually responsible for the evolution of proto¬ 
planetary discs. Our results point toward a basic paradigm of 
protoplanetary disc evolution that includes both laminar and 
turbulent accretion, depending upon the strength and orien¬ 
tation of the magnetic field. The outer disc will be accreting, 
independent of the orientation of the large-scale vertical mag- 
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netic field, through either turbulent or laminar stresses. The 
details of the accretion dynamics depend upon the strength 
of the net vertical magnetic field, and may be testable with 
current sub-millimeter observations ( [Simon et al.||2015l ). On 
smaller radial scales, large-scale correlations in the radial and 
toroidal components of the magnetic field and magnetocen¬ 
trifugal winds l [Bai & Stone||2013a|b[ l are likely to be dom¬ 
inant, although we have emphasized that even here there 
is the potential for turbulence. Possible observational tests 
on ~au scales include the detection of an MHD wind or 
of bi-modality in the nature and amplitude of the angular- 
momentum transport that is predicted to arise from the Hall 
effect. 

In summary, while the nature of accretion in proto¬ 
planetary discs is intimately tied to the details of ioniza¬ 
tion/recombination physics, thereby making precise numer¬ 
ical predictions difficult, there remain robust, qualitative pre¬ 
dictions related to the importance of the large-scale vertical 
magnetic field in driving accretion and the bi-modality of disc 
properties that results from the interaction of the Hall effect 
and this large-scale field. With powerful new observing facil¬ 
ities, such as ALMA, testing these predictions is now entering 
the realm of feasibility. 
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